// eee.Sheffield.PZ.Imaging
//
// Copyright ? Ping Zou, 2007
// sg71.cherub@gmail.com


using System;
using System.Collections.Generic;
using System.Text;
using eee.Sheffield.PZ.Math;
using System.Drawing;
using System.Drawing.Imaging;
using System.IO;

namespace eee.Sheffield.PZ.Imaging.ImageFilter
{
    /// <summary>
    /// GVF return type
    /// </summary>
    public enum GVFReturnType { GVFMagnitude, UMagnitude, VMagnitude }

    public class GradientVectorFlow : ImageFilter//, IGraylevelImageFilter
    {
        #region Fields
        /// <summary>
        /// u part, unit vector
        /// </summary>
        private PZMath_matrix _u = new PZMath_matrix();        
        /// <summary>
        /// v part, unit vector
        /// </summary>
        private PZMath_matrix _v = new PZMath_matrix();        
        /// <summary>
        /// vector length
        /// </summary>
        private PZMath_matrix _vectorMagnitude = new PZMath_matrix();
        /// <summary>
        /// standard deviation of the distribution
        /// </summary>
        private double _sigma = 1.0;
        /// <summary>
        /// regularization parameter mu
        /// </summary>
        private double _mu = 0.2;
        /// <summary>
        /// iteration step, deltat
        /// </summary>
        private double _deltaT;
        /// <summary>
        /// coefficient r
        /// </summary>
        private double _r;   
        private GVFReturnType _GVFReturnType = GVFReturnType.GVFMagnitude;
        //private string _flag = "[Image -> Gradient -> Gradient Vector Flow]";
        #endregion

        #region Properties
        public PZMath_matrix U { get { return _u; } set { _u = value; } }
        public PZMath_matrix V { get { return _v; } set { _v = value; } }
        public PZMath_matrix VectorMagnitude { get { return _vectorMagnitude; } set { _vectorMagnitude = value; } }
        public double Sigma { get { return _sigma; } set { _sigma = value; } }
        public double Mu { get { return _mu; } set { _mu = value; } }
        public double DeltaT { get { return _deltaT; } set { _deltaT = value; } }
        public double R { get { return _r; } set { _r = value; } }
        public GVFReturnType GVFReturnType { get { return _GVFReturnType; } set { _GVFReturnType = value; } }
        //public string Flag() { return _flag; }
        //public string Log() { return _log; }
        #endregion

        #region constructor
        /// <summary>
        /// empty constructor
        /// </summary>
        public GradientVectorFlow() : this(0, 1) { }
        /// <summary>
        /// imput standard deviation of the distribution
        /// </summary>
        /// <param name="sigma"></param>
        public GradientVectorFlow(double sigma, double mu)
        {
            _flag = "[Image -> Gradient -> Gradient Vector Flow]";
            _normalize = true;
            _sigma = sigma;
            _mu = mu;
            _deltaT = 1.0 / 4.0 / _mu / 10;  // 1/2 upper bound (1 / 4 / mu)
            _r = _mu * _deltaT;
        }
        #endregion

        #region Apply Filter
        //public Bitmap ApplyGraylevelImageFilter(Bitmap srcImage) { return Apply(srcImage, out _log); }
        //public PZMath_matrix ApplyDoubleGraylevelImageFilter(PZMath_matrix srcMatrix) { return ApplyDouble(srcMatrix, out _log); }

        //private Bitmap Apply(Bitmap srcImage, out string logText)
        //{
        //    PZMath_matrix srcMatrix = _converter.ConvertGraylevelBitmapToMatrix(srcImage);
        //    PZMath_matrix dstMatrix = ApplyDouble(srcMatrix, out logText);
        //    Bitmap dstImage = _converter.ConvertMatrixToGraylevelBitmap(dstMatrix, true);
        //    return dstImage;
        //}     

        protected override PZMath_matrix ApplyDoubleGraylevelImageFilter(PZMath_matrix srcMatrix, out string logText)
        {
            #region source and dst matrix
            // get source
            int width = srcMatrix.ColumnCount;
            int height = srcMatrix.RowCount;

            // create dstMatrix
            PZMath_matrix dstMatrix = new PZMath_matrix(height, width);
            int heightM1 = height - 1;
            int widthM1 = width - 1;
            int heightM2 = height - 2;
            int widthM2 = width - 2;
            #endregion

            // 1. edge map
            
            EdgeMap edgeMapFilter = new EdgeMap(_sigma);
            PZMath_matrix edgeMap = edgeMapFilter.ApplyDoubleGraylevelImageFilter(srcMatrix);
            //edgeMap.Normalize(0, 1);
            //PZMath_matrix edgeMap = new PZMath_matrix(height, width);
            //for (int y = 0; y < height; y++)
            //    for (int x = 0; x < width; x++)
            //        edgeMap[y, x] = 1 - srcMatrix[y, x] / 255.0;

            //edgeMap.WriteFile("edgemap.txt");

            // 2. gradient x and gradient y of edge map
            edgeMap.Normalize();
            edgeMap = edgeMap.BoundMirrorExpand(3, 3);
            //edgeMap.WriteFile("edgemap.txt");

            int expandedHeight = edgeMap.RowCount;
            int expandedWidth = edgeMap.ColumnCount;

            // 3. gradient of f            
            GradientXCentralDifference gradientXFilter = new GradientXCentralDifference();
            GradientYCentralDifference gradientYFilter = new GradientYCentralDifference();
            PZMath_matrix gradientX = gradientXFilter.ApplyDoubleGraylevelImageFilter(edgeMap);
            PZMath_matrix gradientY = gradientYFilter.ApplyDoubleGraylevelImageFilter(edgeMap);
            string gradientXLog = gradientXFilter.LogText;
            string gradientYLog = gradientYFilter.LogText;
            //gradientX.WriteFile("fx.txt");
            //gradientY.WriteFile("fy.txt");

            // 4. init u, and v to gradientX and gradientY, respectively
            PZMath_matrix U = new PZMath_matrix(expandedHeight, expandedWidth);
            PZMath_matrix V = new PZMath_matrix(expandedHeight, expandedWidth);
            U.MemCopyFrom(gradientX);
            V.MemCopyFrom(gradientY);

            PZMath_matrix sqrMagf = new PZMath_matrix(expandedHeight, expandedWidth);
            for (int y = 0; y < expandedHeight; y++)
                for (int x = 0; x < expandedWidth; x++)
                    sqrMagf[y, x] = gradientX[y, x] * gradientX[y, x] + gradientY[y, x] * gradientY[y, x];
            //sqrMagf.WriteFile("sqrmagf.txt");

            // 5. iteration
            PZMath_matrix laplacianU = new PZMath_matrix(expandedHeight, expandedWidth);
            PZMath_matrix laplacianV = new PZMath_matrix(expandedHeight, expandedWidth);
            DiscreteLaplacian laplacianFilter = new DiscreteLaplacian();
            int totalIteration = (int)System.Math.Sqrt(System.Math.Max(srcMatrix.RowCount, srcMatrix.ColumnCount));
            totalIteration = 100;
            for (int iteration = 0; iteration < totalIteration; iteration++)
            {
                U.BoundMirrorEnsure(3, 3);
                V.BoundMirrorEnsure(3, 3);
                //U.WriteFile("u.txt");
                //V.WriteFile("v.txt");
                
                laplacianU = laplacianFilter.ApplyDoubleGraylevelImageFilter(U);
                laplacianV = laplacianFilter.ApplyDoubleGraylevelImageFilter(V);
                //laplacianU.WriteFile("laplacianu.txt");
                //laplacianV.WriteFile("laplacianv.txt");

                for (int y = 0; y < expandedHeight; y++)
                {
                    for (int x = 0; x < expandedWidth; x++)
                    {
                        U[y, x] = U[y, x] + _mu * 4.0 * laplacianU[y, x] - sqrMagf[y, x] * (U[y, x] - gradientX[y, x]);
                        V[y, x] = V[y, x] + _mu * 4.0 * laplacianV[y, x] - sqrMagf[y, x] * (V[y, x] - gradientY[y, x]);
                    }
                }
                //U.WriteFile("u.txt");
                //V.WriteFile("v.txt");
            }

            // 6. boundary shrinnk
            _u = U.BoundMirrorShrink(3, 3);
            _v = V.BoundMirrorShrink(3, 3);

            // 7. vector length and unit vector
            _vectorMagnitude = new PZMath_matrix(height, width);
            for (int y = 0; y < height; y++)
            {
                for (int x = 0; x < width; x++)
                {
                    _vectorMagnitude[y, x] = System.Math.Sqrt(_u[y, x] * _u[y, x] + _v[y, x] * _v[y, x]);
                    _u[y, x] = _u[y, x] / (_vectorMagnitude[y, x] + PZMath_machine.PZMath_DBL_MIN);
                    _v[y, x] = _v[y, x] / (_vectorMagnitude[y, x] + PZMath_machine.PZMath_DBL_MIN);
                }
            }

            logText = "succeed.";
            // return GVF images
            if (_GVFReturnType == GVFReturnType.GVFMagnitude)
            {
                _vectorMagnitude.Normalize(0.0, 1.0);
                return _vectorMagnitude;
            }
            else if (_GVFReturnType == GVFReturnType.VMagnitude)
            {
                _v.Normalize(0.0, 1.0);
                return _v;
            }
            else
            {
                _u.Normalize(0.0, 1.0);
                return _u;
            }
        } // ApplyDoubleGraylevelImageFilter()

        protected override PZMath_matrix[] ApplyDoubleColourImageFilter(PZMath_matrix[] srcMatrixArray, out string logText)
        {
            // TODO: add this function

            throw new Exception("The method or operation is not implemented.");
        } // ApplyDoubleColourImageFilter()
        #endregion

        #region I/O
        /// <summary>
        /// Tecplot IJ-ordered ascii data file
        /// </summary>
        /// <param name="fileName"></param>
        public void WriteGVFTecplotFile(string fileName)
        {
            FileStream fileStream = new FileStream(fileName, FileMode.Create, FileAccess.Write);
            StreamWriter writer = new StreamWriter(fileStream);
            int height = _u.RowCount;
            int width = _v.ColumnCount;
            writer.WriteLine("TITLE     = " + fileName);
            writer.WriteLine("VARIABLES = \"x\"");
            writer.WriteLine("\"y\"");
            writer.WriteLine("\"u\"");
            writer.WriteLine("\"v\"");
            writer.WriteLine("\"magenitude\"");
            writer.WriteLine("ZONE T = \"Zone 1\"");
            writer.WriteLine("I = " + width + ", J=" + height + ", ZONETYPE=Ordered");
            writer.WriteLine("DATAPACKING = POINT");
            for (int y = 0; y < height; y++)
            {
                for (int x = 0; x < width; x++)
                {
                    writer.WriteLine(String.Format("{0,-20:0.00}{1,-20:0.00}{2,-20:0.00}{3,-20:0.00}{4,-20:0.00}",
                        x, y, _u[y, x], _v[y, x], _vectorMagnitude[y, x]));
                }
            }
            writer.Flush();
            writer.Close();
            fileStream.Close();
        }

        public void WriteGVFonCentreLineTecplotFile(string fileName, Bitmap centerLineImage)
        {
            WriteGVFonCentreLineTecplotFile(fileName, centerLineImage, 0);

            //FileStream fileStream = new FileStream(fileName, FileMode.Create, FileAccess.Write);
            //StreamWriter writer = new StreamWriter(fileStream);
            //int height = _u.RowCount;
            //int width = _v.ColumnCount;
            //writer.WriteLine("TITLE     = " + fileName);
            //writer.WriteLine("VARIABLES = \"x\"");
            //writer.WriteLine("\"y\"");
            //writer.WriteLine("\"u\"");
            //writer.WriteLine("\"v\"");
            //writer.WriteLine("\"magenitude\"");
            //writer.WriteLine("ZONE T = \"Zone 1\"");
            //writer.WriteLine("I = " + width + ", J=" + height + ", ZONETYPE=Ordered");
            //writer.WriteLine("DATAPACKING = POINT");
            //for (int y = 0; y < height; y++)
            //{
            //    for (int x = 0; x < width; x++)
            //    {
            //        if (Convert.ToInt32(centerLineImage.GetPixel(x, y).R) == 255)
            //        // no center line
            //        {
            //            writer.WriteLine(String.Format("{0,-20:0.00}{1,-20:0.00}{2,-20:0.00}{3,-20:0.00}{4,-20:0.00}",
            //                x, y, 0, 0, 0));
            //        }
            //        else
            //        {
            //            writer.WriteLine(String.Format("{0,-20:0.00}{1,-20:0.00}{2,-20:0.00}{3,-20:0.00}{4,-20:0.00}",
            //                x, y, _u[y, x], _v[y, x], _vectorMagnitude[y, x]));
            //        }
            //    }
            //}
            //writer.Flush();
            //writer.Close();
            //fileStream.Close();
        }

        public void WriteGVFonCentreLineTecplotFile(string fileName, Bitmap centerLineImage, int offset)
        {
            FileStream fileStream = new FileStream(fileName, FileMode.Create, FileAccess.Write);
            StreamWriter writer = new StreamWriter(fileStream);
            int height = _u.RowCount;
            int width = _v.ColumnCount;
            writer.WriteLine("TITLE     = " + fileName);
            writer.WriteLine("VARIABLES = \"x\"");
            writer.WriteLine("\"y\"");
            writer.WriteLine("\"u\"");
            writer.WriteLine("\"v\"");
            writer.WriteLine("\"magenitude\"");
            writer.WriteLine("ZONE T = \"Zone 1\"");
            writer.WriteLine("I = " + width + ", J=" + height + ", ZONETYPE=Ordered");
            writer.WriteLine("DATAPACKING = POINT");
            for (int y = 0; y < height; y++)
            {
                for (int x = 0; x < width; x++)
                {
                    if (Convert.ToInt32(centerLineImage.GetPixel(x, y).R) == 255)
                    // no center line
                    {
                        writer.WriteLine(String.Format("{0,-20:0.00}{1,-20:0.00}{2,-20:0.00}{3,-20:0.00}{4,-20:0.00}",
                            x, y, 0, 0, 0));
                    }
                    else
                        // draw offset line
                    {                        
                        writer.WriteLine(String.Format("{0,-20:0.00}{1,-20:0.00}{2,-20:0.00}{3,-20:0.00}{4,-20:0.00}",
                            x, y, _u[y + offset, x], _v[y + offset, x], _vectorMagnitude[y + offset, x]));
                    }
                }
            }
            writer.Flush();
            writer.Close();
            fileStream.Close();
        }

        public void WriteFileU(string fileName)
        {
            _u.WriteFile(fileName);
        }

        public void WriteFileV(string fileName)
        {
            _v.WriteFile(fileName);
        }

        public void WriteFileVectorMagnitude(string fileName)
        {
            _vectorMagnitude.WriteFile(fileName);
        }

        public void ReadFileU(string fileName)
        {
            _u.ReadFile(fileName);
        }

        public void ReadFileV(string fileName)
        {
            _v.ReadFile(fileName);
        }

        public void ReadFileVectorMagnitude(string fileName)
        {
            _vectorMagnitude.ReadFile(fileName);
        }
        #endregion
    }
}